Gain Mismatch systematics

In [23]:
import healpy as hp 
import pylab as pl 


mapdic={0:'T',1:'Q', 2:'U'}

Goal of these studies

We study the systematics effects induced by gain mismatches for Litebird scanning strategy and a fake focal plane with just 2 detector pairs (to save computation time and memory at NERSC).

We inject gain drifts only in one of the two detectors as random Gaussian fluctations around 0 and with a 5% width during 1 precession period. In our last In our last post, we've described the adopted methodology to estimate gain drifts with a template regression.

Simulation Settings

  • simulate 1 yr of observations with LB scanning strategy parameters, $\alpha= 45$, $\beta =50$, $spinrate=0.1 rpm $, $precperiod=96.5 min$
  • Focal plane encodes 4 detectors (two det. pairs) at 140 GHz , FWHM= 25 arcmin, NET= 38.6 $\mu K \sqrt{s}$, sample rate= 40 Hz
  • not spinning HWP /!\
  • $f_{knee}$ is allowed to range 0,20, 50 mHz
  • Input signal : CMB T only map + Dipole + Noise
  • Output maps : destriped with MADAM (baseline=1min, noisefilter=Off)

Choosing the template

In this posting we firstly want to quantify how the gain correction procedure is affected by the template chosen to perform the regression . We thus consider signal only simulations encoding as a template for the regression the CMB T only maps w/ and w/o the Orbital Dipole.

We allow the gain to drift during 1 precession period onlyin one of the two detectors in each pair ( namely Gaussian fluctations around 0 and with a 5% width).

We estimate the T->P leakage in terms of power spectra of residual maps, $ \Delta m $ , defined as

$ \Delta m = m_{out} - m_{in} $

In [30]:
cmbmaps=hp.read_map(field=range(3), filename='/Users/peppe/work/satellite_sims/sigonly/CMB_TQU_256_r0_1deg.fits',
                    nest=True,verbose =False)

lmax=3*hp.get_nside(cmbmaps[0])
l=pl.arange(2, lmax+1)
bl=hp.gauss_beam(lmax=lmax, fwhm=pl.deg2rad(10/60))[2:]
cmbmaps=hp.read_map(field=range(3), filename='/Users/peppe/work/satellite_sims/sigonly/CMB_TQU_256_r0_1deg.fits',
                    nest=True,verbose =False)

inputclfile ='/Users/peppe/work/satellite_sims/sigonly/inputcl.fits' 
try :
    cl = hp.read_cl(filename=inputclfile)
except IOError:
    cl=hp.anafast(cmbmaps, lmax=lmax)
    hp.write_cl(cl=cl, filename='inputclfile' ) 
hitmap=hp.read_map('/Users/peppe/work/satellite_sims/madam_hmap.fits',verbose=False )
hp.mollview(pl.log10( hitmap), min=3.5, max=4, unit='Log (hitcounts)', sub=111 , title='' )

Regressing w/ Dipole

Let's have a look at the input maps

In [47]:
dipolemap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/dipole/madam_bmap.fits', 
                       verbose=False)

sigonlymap = hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/dipole_CMB_T/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(5,5))
hp.mollview(sigonlymap [0], unit='K ', title='Dipole T  ', sub=221 , fig=fig )
cmbT =sigonlymap [0] - dipolemap[0]
hp.mollview(cmbT , unit='K ', title='CMB T  ', sub=222 , fig=fig )
hp.mollview(sigonlymap [1], unit='K ', title='Q  ', sub=223 , fig=fig )
hp.mollview(sigonlymap [2], unit='K ', title='U  ', sub=224 , fig=fig )

As expected the Dipole amplitude is 3 mK $\sim 10\, $times larger than CMB Temperature.

NB: Polariz. maps are essentially null.

Injecting gain drifts

We look at residual maps

In [48]:
driftmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(  driftmap [0] -sigonlymap[0] , unit='K ', title=r'$\Delta T$  ', sub=131 , fig=fig   )
hp.mollview(driftmap [1], unit='K ', title=r'$\Delta Q$', sub=132 , fig=fig )
hp.mollview(driftmap [2], unit='K ', title=r'$\Delta U$', sub=133 , fig=fig  )
cldrifts=hp.anafast(driftmap-sigonlymap, lmax=lmax)

We get non null signal on Q, and U maps: that's the expected T-> P leakage from gain drifts. Moreover, all the residuals are at the very same order of magnitude.
A ring-like pattern is observed in all maps and seems to fall exactly where the dipole changes sign. To check this, we select pixels in the dipole maps whose value is in the interval $\left[ -800, 800 \right] \mu K$.

In [54]:
ring=pl.ma.masked_inside(dipolemap[0], -8e-4, 8e-4) 
mask=pl.zeros_like(dipolemap[0])  
mask[ring.mask]=1. 
fig=pl.figure(figsize=(15,8))
hp.mollview(mask*( driftmap [0] -sigonlymap[0] ),  min=-3e-6, max=3e-6 ,
              sub=131, title=' T')

hp.mollview(mask*( driftmap [1] -sigonlymap[1] ),  min=-3e-6, max=3e-6 ,
              sub=132, title=' Q')
hp.mollview(mask*( driftmap [2] -sigonlymap[2] ),   min=-3e-6, max=3e-6 ,
              sub=133,title='  U')
 

We plot above the orthographic view centered in the north Ecliptic pole of the residual maps we 've seen above. We masked out all the residual pixels where the Dipole intensity is outside the selection interval and we notice that the ring pattern perfectly coincides to the one defined from the dipole (the circles are parallel ).

We can therefore conclude that the ring feature is introduced by the dipole

Correcting for gain drifts

We now estimate the gains by regressing against the Dipole template.

The correction is applied at the TOD level.

we run map-making on top of the post-correction TODs.

In [34]:
driftcorrmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/corr/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview((driftcorrmap [0]  -sigonlymap[0])    ,unit='K', title=r'$\Delta T  $  ', sub=131 , 
            fig=fig  )
hp.mollview(driftcorrmap [1] , unit='K',title=r'$\Delta Q  $', sub=132 , fig=fig ,
            )
hp.mollview(driftcorrmap [2] ,unit='K', title=r'$\Delta U   $', sub=133 , fig=fig 
             )
cldriftcorr=hp.anafast(driftcorrmap-sigonlymap , lmax=lmax)
In [35]:
print "Pre-Correction\nMap \t Mean \t \t Std "
for i  in xrange(3): 
    print mapdic[i], (driftmap[i] - sigonlymap[i]).mean()  ,(driftmap[i] - sigonlymap[i]).std()  
print "\nPost-correction "
for i  in xrange(3): 
    print mapdic[i], (driftcorrmap[i] - sigonlymap[i]).mean()  ,(driftcorrmap[i] - sigonlymap[i]).std()  
Pre-Correction
Map 	 Mean 	 	 Std 
T 1.5814402163034245e-07 4.0165867310115265e-06
Q -1.3173931719076805e-09 5.21026343757933e-06
U 1.342745920389484e-08 5.210053553479623e-06

Post-correction 
T 1.7448183140128026e-08 1.9795333558065137e-06
Q -4.902921210001606e-10 3.7365168916884247e-07
U -1.4399108583470596e-09 3.629386507584052e-07

You may notice that residual Q and U maps are a factor ~10 times smaller than the ones w/o any gain correction (see above ). Together with a dipole-like residual global structure, we observe a couple of scanning strategy features (the two horizontal stripes along the Ecliptic plane and in the Ecliptic poles) . These are the regions where pixels are observed $\sim 3 $ times more wrt the rest of the sky, as shown in the hitmap above.

We will notice below that having residuals that resemble the template adopted to perform the regression is a feature related to this methodology. By dividing the residual maps by the dipole map $D$, we get the fractional residual maps and they shouldn't present any feature similar to the template.

In [220]:
fig=pl.figure(figsize=(15,8))
hp.mollview( (driftcorrmap [0]  -sigonlymap[0])/ (sigonlymap[0])   , title=r'$\Delta T /D$  ', sub=131 , 
            fig=fig , min=-0.01,max=0.01 )
hp.mollview(driftcorrmap [1]/ sigonlymap[0] , title=r'$\Delta Q /D$', sub=132 , fig=fig ,
            min=-3e-4,max=3e-4 )
hp.mollview(driftcorrmap [2]/ sigonlymap[0] , title=r'$\Delta U / D $', sub=133 , fig=fig 
            ,min=-3e-4,max=3e-4 )

As expected, the fractional residuals present no more dipole-like features (the ring of pixels is now narrower and is due to $D$ being close to zero in those pixels. We set the colorbar in such a way that fractional differences in polarization maps are at most $0.03 \%$ smaller than the dipole, meaning a residual map fluctuating $\lesssim 1 \mu K$. For the fractional residual temperature map we set as a reference residual value the $1 \%$ of the dipole, i.e. $\sim 30 \mu K $.

We then compute full sky power spectra (with anafast) from the residual maps, $\Delta m$.

In [189]:
fig=pl.figure(figsize=(16, 5))
pl.subplot(131)
pl.title('TT')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[0][2:]/bl**2 *1e12,'k', label='CMB')
pl.loglog(l, cldrifts[0][2:]/bl**2 *1e12, label='gain drifts w/o noise ')
pl.loglog(l, cldriftcorr[0][2:]/bl**2 *1e12, label='gain correct. w/o noise ')
pl.legend(loc='upper right')
pl.xlim([2,600])

pl.subplot(132)
pl.title('EE')
pl.loglog(l,cl[1][2:]/bl**2  *1e12,'k')
pl.loglog(l, cldrifts[1][2:]/bl**2 *1e12)
pl.loglog(l, cldriftcorr[1][2:]/bl**2 *1e12)

pl.ylim([1e-8,1e0])
pl.xlim([2,600])
pl.subplot(133)
pl.title('BB')
pl.loglog(l, cl[2][2:] /bl**2 *1e12,'k')
pl.loglog(l, cldrifts[2][2:]/bl**2 *1e12)
pl.loglog(l, cldriftcorr[2][2:]/bl**2 *1e12)
pl.loglog(l, pl.ones_like(l)*4.87e-6,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e0])
foo=pl.xlim([2,600])

Solid blue (orange) are power spectra of pre(post)-gain-correction residual maps. Since no noise has been coadded to the TODs so far, the plots above represent the level of accuracy our estimator might achieve when the regression exploits the dipole as a template. As a reference, we plot the CMB spectra from $\Lambda CDM$. We also include in BB plot the solid red line resembling the ultimate Litebird whitenoise level of $7.6 \mu $K arcmin, after FG subtraction.

Pros: By comparing the blue with the orange solid lines should highlights the effects of gain correction. As expected they are small for the temperature maps, but in terms of EE and BB power the correction is very effective: it may correct the T->P leakage as much as 3 orders of magnitude in the $\sim$ degree scales.

Cons: the correction doesn't work at $\ell<10$ scales.

Regressing w/o Dipole

we keep the same configuration as before, we just changed the input maps encoding now only CMB T.

In [36]:
sigonlymap = hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/CMB_T/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(12,8))
hp.mollview(sigonlymap[0]  , unit='K ', title=' T  ', sub=131 , fig=fig )
hp.mollview(sigonlymap [1] , unit='K ', title='Q  ', sub=132 , fig=fig  )
hp.mollview(sigonlymap [2], unit='K ', title='U  ', sub=133, fig=fig )
In [37]:
driftmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/wodipole/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(  driftmap [0] -sigonlymap[0] , unit='K ', title=r'$\Delta T$  ', sub=131 , fig=fig   )
hp.mollview(driftmap [1], unit='K ', title=r'$\Delta Q$', sub=132 , fig=fig )
hp.mollview(driftmap [2], unit='K ', title=r'$\Delta U$', sub=133 , fig=fig  )
cldrifts=hp.anafast(driftmap-sigonlymap, lmax=lmax)
                       

The pre-correction residual maps look as noise fluctuating around $10^{-6}$. This also confirms that the ring-like residual is introduced by the dipole.

We now apply the gain correction and output the maps.

In [38]:
driftcorrmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/corr/wodipole/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(driftcorrmap [0]- sigonlymap[0] , unit='K ', title=r'$\Delta T$  ', sub=131 , fig=fig, norm='hist' )
hp.mollview(driftcorrmap [1] , unit='K ', title=r'$\Delta Q$', sub=132 , fig=fig, norm='hist' )
hp.mollview(driftcorrmap [2], unit='K ', title=r'$\Delta U$', sub=133 , fig=fig, norm='hist' )
cldriftcorr=hp.anafast(driftcorrmap-sigonlymap , lmax=lmax)
In [39]:
print "Pre-Correction\nMap \t Mean \t \t Std "
for i  in xrange(3): 
    print mapdic[i], (driftmap[i] - sigonlymap[i]).mean()  ,(driftmap[i] - sigonlymap[i]).std()  
print "\nPost-correction "
for i  in xrange(3): 
    print mapdic[i], (driftcorrmap[i] - sigonlymap[i]).mean()  ,(driftcorrmap[i] - sigonlymap[i]).std()  
Pre-Correction
Map 	 Mean 	 	 Std 
T -1.3031631041792773e-09 1.4703558705737005e-07
Q -1.2506098329267013e-10 2.0635097291731076e-07
U 2.076359928571984e-09 2.0687227199345813e-07

Post-correction 
T 4.679025172596079e-09 4.970213410028304e-08
Q -1.4803030274676635e-09 1.6309438817575164e-08
U -1.4659737498610382e-09 1.5764294530864922e-08

Post-correction residual maps are one order of magnitude smaller than the pre-correction ones, indicating that the methodology is working very well. Also, we achieve 10 times lower residual by regressing with CMB T maps than in the Dipole case.

Moreover, note that to highlight differences we compress the colorscale. A pattern proportional to a Temperature map is also visible in all the residual maps, as stated above as well in case of Dipole.
Scanning strategy feature are also visible in Q and U residual maps. We thus divide the the residual maps by the input CMB temperature one to assess fractional residuals.

In [40]:
fig=pl.figure(figsize=(15,8))
hp.mollview( (driftcorrmap [0]  -sigonlymap[0])/sigonlymap[0]    , title=r'$\Delta T /T$  ', sub=131 , 
            fig=fig  , norm='hist')
hp.mollview(driftcorrmap [1]/ sigonlymap[0] , title=r'$\Delta Q /T$', sub=132 , fig=fig ,
            min=-1e-3,max=1e-3 )
hp.mollview(driftcorrmap [2]/ sigonlymap[0] , title=r'$\Delta U /T $', sub=133 , fig=fig 
            ,min=-1e-3,max=1e-3 )

As already stated above, the Q and U fractional residuals presents just the scanning strategy features. In particular, these residuals are below the $0.1 \%$ level wrt to T anisotropies ($\sim 0.1 \mu K$) .

For Temperature fractional residuals they all fluctuates on a positive value ranging from $0.04 \%$ to $0.1 \%$

In [226]:
fig=pl.figure(figsize=(16, 5))
pl.subplot(131)
pl.title('TT')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[0][2:]/bl**2 *1e12,'k', label='CMB')
pl.loglog(l, cldrifts[0][2:]/bl**2 *1e12, label='gain drifts w/o noise ')
pl.loglog(l, cldriftcorr[0][2:]/bl**2 *1e12, label='gain correct. w/o noise ')
pl.legend(loc='upper right')
pl.xlim([2,600])

pl.subplot(132)
pl.title('EE')
pl.loglog(l,cl[1][2:]/bl**2  *1e12,'k')
pl.loglog(l, cldrifts[1][2:]/bl**2 *1e12)
pl.loglog(l, cldriftcorr[1][2:]/bl**2 *1e12)

pl.ylim([1e-8,1e0])
pl.xlim([2,600])
pl.subplot(133)
pl.title('BB')
pl.loglog(l, cl[2][2:] /bl**2 *1e12,'k')
pl.loglog(l, cldrifts[2][2:]/bl**2 *1e12)
pl.loglog(l, cldriftcorr[2][2:]/bl**2 *1e12)
pl.loglog(l, pl.ones_like(l)*4.87e-6,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e0])
foo=pl.xlim([2,600])

The level of systematics introduced by gain drifts in the case of T only sims, is lower compared to the one with T+Dipole, though even large scale B-modes ($\ell <10 $) of post-correction maps are as much as comparable as the expected level of lensing B-modes at those angular scales. However, we are safely below the white noise level after Foreground subtraction (solid red line).

Injecting faster (10 min) gain drifts

So far, we've studied the case of a gain slowly drifting within a time scale of $\sim 1.5 $ h. In order to be more realistic we considered a shorter time scale, i.e. $10$ min.

In [42]:
driftmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/wodipole/10min/8det/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(  driftmap [0] -sigonlymap[0] , unit='K ', title=r'$\Delta T$  ', sub=131 , fig=fig  )
hp.mollview(driftmap [1], unit='K ', title=r'$\Delta Q$', sub=132 , fig=fig   )
hp.mollview(driftmap [2], unit='K ', title=r'$\Delta U$', sub=133 , fig=fig  )
cldrifts=hp.anafast(driftmap-sigonlymap, lmax=lmax)
                       

We get similar residual amplitudes as in the slowly varying gain drift case. Gain drifts injected into simulations behave as an extra-white noise to the maps (In facts, at the power spectrum level the pre-correction power spectra (solid blue lines) are almost flat ).

We now apply the gain correction and output the maps.

In [43]:
driftcorrmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/corr/wodipole/10min/8det/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(driftcorrmap [0]- sigonlymap[0] , unit='K ', title=r'$\Delta T$  ', sub=131 , fig=fig, norm='hist' )
hp.mollview(driftcorrmap [1] , unit='K ', title=r'$\Delta Q$', sub=132 , fig=fig, norm='hist' )
hp.mollview(driftcorrmap [2], unit='K ', title=r'$\Delta U$', sub=133 , fig=fig, norm='hist' )
cldriftcorr=hp.anafast(driftcorrmap-sigonlymap , lmax=lmax)
In [44]:
print "Pre-Correction\nMap \t Mean \t \t Std "
for i  in xrange(3): 
    print mapdic[i], (driftmap[i] - sigonlymap[i]).mean()  ,(driftmap[i] - sigonlymap[i]).std()  
print "\nPost-correction "
for i  in xrange(3): 
    print mapdic[i], (driftcorrmap[i] - sigonlymap[i]).mean()  ,(driftcorrmap[i] - sigonlymap[i]).std()  
Pre-Correction
Map 	 Mean 	 	 Std 
T -4.13637445375782e-10 1.0350415623196197e-07
Q 3.025793933049551e-10 1.4599027084427707e-07
U 9.388814206877955e-11 1.4670392014577685e-07

Post-correction 
T 4.876082384539356e-09 4.9387105564556566e-08
Q -1.4963132550693472e-09 1.520816262285405e-08
U -1.4481320061591677e-09 1.4783103501021382e-08
In [45]:
fig=pl.figure(figsize=(16, 5))
pl.subplot(131)
pl.title('TT')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[0][2:]/bl**2 *1e12,'k', label='CMB')
pl.loglog(l, cldrifts[0][2:]/bl**2 *1e12, label='gain drifts w/o noise ')
pl.loglog(l, cldriftcorr[0][2:]/bl**2 *1e12, label='gain correct. w/o noise ')
pl.legend(loc='upper right')
pl.xlim([2,600])

pl.subplot(132)
pl.title('EE')
pl.loglog(l,cl[1][2:]/bl**2  *1e12,'k')
pl.loglog(l, cldrifts[1][2:]/bl**2 *1e12)
pl.loglog(l, cldriftcorr[1][2:]/bl**2 *1e12)

pl.ylim([1e-8,1e0])
pl.xlim([2,600])
pl.subplot(133)
pl.title('BB')
pl.loglog(l, cl[2][2:] /bl**2 *1e12,'k')
pl.loglog(l, cldrifts[2][2:]/bl**2 *1e12)
pl.loglog(l, cldriftcorr[2][2:]/bl**2 *1e12)
pl.loglog(l, pl.ones_like(l)*4.87e-6,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e0])
foo=pl.xlim([2,600])

Map residuals are comparable to the ones with the slowly varying gain fluctuations. This is somewhat expected as long as we estimate the gain variations with the very same time-scale the variations has been injected.

Regressing with T and E scalar anisotropies (no lensing B-modes)

To get more realistic simulations we include $E$ -modes polarization to the signal to be observed. So that, now pair differences encode a non-null signal. We keep 10 minutes as the gain drift time scale.

In [227]:
sigonlymap = hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/scalar_only/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(14,12))
hp.mollview(sigonlymap[0]  , unit='K ', title=' T  ', sub=131 , fig=fig )
hp.mollview(sigonlymap [1] , unit='K ', title='Q  ', sub=132 , fig=fig   )
hp.mollview(sigonlymap [2], unit='K ', title='U  ', sub=133 , fig=fig )
In [228]:
driftmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/uncorr/scalar_only/10min/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(  driftmap [0] -sigonlymap[0] , unit='K ', title=r'$\Delta T$  ', sub=131 , fig=fig , norm='hist' )
hp.mollview(driftmap [1]-sigonlymap[1], unit='K ', title=r'$\Delta Q$', sub=132 , fig=fig, norm='hist')
hp.mollview(driftmap [2] -sigonlymap[2], unit='K ', title=r'$\Delta U$', sub=133 , fig=fig, norm='hist' )
cldrifts=hp.anafast(driftmap-sigonlymap, lmax=lmax)
In [229]:
driftcorrmap= hp.read_map(field=range(3) , 
                       filename='/Users/peppe/work/satellite_sims/sigonly/corr/scalar_only/10min/madam_bmap.fits', 
                       verbose=False)
fig=pl.figure(figsize=(15,8))
hp.mollview(driftcorrmap [0]- sigonlymap[0] , unit='K ', title=r'$\Delta T$  ', sub=131 , fig=fig, norm='hist' )
hp.mollview(driftcorrmap [1]- sigonlymap[1], unit='K ', title=r'$\Delta Q$', sub=132 , fig=fig, norm='hist' )
hp.mollview(driftcorrmap [2]- sigonlymap[2], unit='K ', title=r'$\Delta U$', sub=133 , fig=fig , norm='hist')
cldriftcorr=hp.anafast(driftcorrmap-sigonlymap , lmax=lmax)
In [46]:
fig=pl.figure(figsize=(16, 5))
pl.subplot(221)
pl.title('TT')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l, cl[0][2:]/bl**2 *1e12,'k', label='CMB')
pl.loglog(l, cldrifts[0][2:]/bl**2 *1e12, label='gain drifts w/o noise ')
pl.loglog(l, cldriftcorr[0][2:]/bl**2 *1e12, label='gain correct. w/o noise ')
pl.legend(loc='upper right')
pl.xlim([2,700])

pl.subplot(222)
pl.title('TE')
s= pl.sign(cl[3][2:])
cycle = pl.rcParams['axes.prop_cycle'].by_key()['color']
pl.loglog(l, cl[3][2:]/bl**2 *1e12,'k' )
pl.loglog(l, -cl[3][2:]/bl**2 *1e12,'k:' )
pl.loglog(l, cldrifts[3][2:]/bl**2 *1e12,color=cycle[0] )
pl.loglog(l, -cldrifts[3][2:]/bl**2 *1e12,'.-',color=cycle[0] )

pl.loglog(l, cldriftcorr[3][2:]/bl**2 *1e12,color=cycle[1] )
pl.loglog(l, -cldriftcorr[3][2:]/bl**2 *1e12,'.',color=cycle[1] )

pl.xlim([2,700])

pl.subplot(223)
pl.title('EE')
pl.ylabel(r'$C_{\ell} [ \mu K^2 ]$', fontsize=16)
pl.loglog(l,cl[1][2:]/bl**2  *1e12,'k')
pl.loglog(l, cldrifts[1][2:]/bl**2 *1e12)
pl.loglog(l, cldriftcorr[1][2:]/bl**2 *1e12)

pl.ylim([1e-8,1e0])
pl.xlim([2,600])
pl.subplot(224)
pl.title('BB')
pl.loglog(l, cl[2][2:] /bl**2 *1e12,'k')
pl.loglog(l, cldrifts[2][2:]/bl**2 *1e12)
pl.loglog(l, cldriftcorr[2][2:]/bl**2 *1e12)
pl.loglog(l, pl.ones_like(l)*4.87e-6,'r', alpha=.5, linewidth=3)
pl.ylim([1e-8,1e-4])
foo=pl.xlim([2,600])

Having now non null E-mode polarization, yields to a non null TE correlation. We therefore show TE spectra with solid (dotted) indicating positive (negative) values. The residual (post-correction) distortion in TE are about 3 orders of magnitude smaller than the TE CMB sepctrum at all multipoles. Since we expect the gain drifts leaking from T to P (no E->B leakage) , the residuals in E and B mode spectra are as comparable as the case presented above encoding just T CMB maps.